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ABSTRACT 

We present the DIRTY radiative transfer model in this paper and a companion paper. This model 
computes the polarized radiative transfer of photons from arbitrary distributions of stars through arbi- 
trary distributions of dust using Monte Carlo techniques. The dust re-emission is done self-consistently 
with the dust absorption and scattering and includes all three important emission paths: equilibrium 
thermal emission, non-equilibrium thermal emission, and the aromatic features emission. The algorithm 
used for the radiative transfer allows for the efficient computation of the appearance of a model system 
as seen from any viewing direction. We present a simple method for computing an upper limit on the 
output quantity uncertainties for Monte Carlo radiative transfer models which use the weighted photon 
approach. 

Subject headings: radiative transfer 



1. INTRODUCTION 

Dust radiative transfer models are crucial to interpret- 
ing many astronomical observations ranging from reflec- 
tion nebulae with a single star to galaxies comprised of 
stars, gas, and dust arranged in complex geometries. Al- 
most all observations are affected by the presence of dust 
intrinsic to the system under study or present in the fore- 
ground. For dust present in the foreground as is the case 
for observations of single stars, the effects of dust can be 
accurately taken into account using a screen geometry. In 
this geometry, the effects are well characterized by an ex- 
tinction curve which is determined by the amount and in- 
trinsic properties (e.g., sizes and materials) of the dust 
(Clayton et al. 2000). However, when the dust is near or 
mixed with the system being studied, its effects become 
dependent not only on the intrinsic properties of the dust 
grains, but on the geometry of the stars and dust in the 
system as well (Witt, Thronson, & Capuano 1992; Witt 
& Gordon 1996, 2000). Since the effects of mixing stars 
and dust are non-trivial, a radiative transfer model is re- 
quired. Using such a model produces a deeper understand- 
ing of the physical properties of the system being studied. 
This applies to a wide variety of astrophysical systems, 
from reflection nebulae to galaxies. There are a variety 
of numerical methods other than Monte Carlo techniques 
which can be used to solve the radiative transfer through 
dust. One example is the doubling method (Zubko & Laor 
2000). While the other methods can achieve faster com- 
putations for specific cases, Monte Carlo techniques are 
the most flexible and allow the study of a wider range of 
astrophysical systems with arbitrary geometries. 

There exist many Monte Carlo dust radiative transfer 
models which have been constructed to study a variety of 
problems. In most models, the computational needs of the 
Monte Carlo technique have been reduced by exploiting 



symmetries and/or making approximations. This neces- 
sarily restricts the variety of systems to which a particu- 
lar model can be successfully applied. Reflection nebulae 
and circumstellar shells have been successfully approxi- 
mated as spherically and azimuthally symmetric systems 
(Witt 1977; Witt et al. 1982; Voshchinnikov & Karjukin 
1994; Calzetti et al. 1995; Witt & Gordon 1996). Models 
for bipolar nebulae, associated with young stellar objects 
and active galactic nuclei, usually exploit the azimuthal 
symmetry of such systems (Yusef-Zadeh, Morris, & White 
1984; Whitney & Hartmann 1992; Fischer, Henning, & 
Yorke 1994; Wood et al. 1998; Young 2000). The gen- 
eral behavior of dust in galaxies has been studied using the 
concept of spherical galactic environments (Witt, Thron- 
son, & Capuano 1992; Witt & Gordon 2000). Sphal 
galaxies have been studied with models using the approx- 
imation of disk galaxies (Bianchi, Ferrara, & Giovanardi 
1996; Wood & Jones 1997; Kuchinski et al. 1998). 

The importance of the multi-phase nature of dust dis- 
tributions has motivated the construction of models which 
can handle arbitrary distributions of dust. Models which 
can handle such dust distributions include those of Wolf, 
Fischer, & Pfau (1998), Varosi & Dwek (1999), Wood 
& Reynolds (1999), Bianchi et al. (2000), and our 
model. More recently, the importance of the re-emission 
of the energy absorbed by dust has been acknowledged 
and increasingly added to some radiative transfer mod- 
els. The dust re-emission can be broken into three com- 
ponents; thermal equilibrium emission (large particles), 
non-equilibrium thermal emission (small particles), and 
the aromatic features (PAH-like particle emission between 
3.3 and ^^20 /im). The inclusion of the dust re-emission to 
radiative transfer models has been done with a range of ap- 
proximations. The easiest approximation is to include dust 
in thermal equilibrium (Varosi & Dwek 1999; Wolf, Hen- 
ning, & Stecklum 1999; Bianchi, Davies, & Alton 2000). 
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Including the non-equilibrium emission and aromatic fea- 
ture emission is much harder and, as a result, requires a 
significant amount of CPU time. Due to the complexity of 
the non-equilibrium and aromatic feature emission, some 
models use a simplified (Silva ct al. 1998) and/or empir- 
ical approach to including them. The model described in 
this paper and a companion paper (Misselt et al. 2000a) 
can handle arbitrary dust distributions and fully includes 
all three dust emission components self-consistently. 

The DIRTY niodc;l (DustI Radiative Transfer, Yeah!) 
grew out of spherically symmetric models which were used 
to study reflection nebulae (Witt et al. 1982; Gordon 
et al. 1994; Calzetti et al. 1995). Using Monte Carlo 
techniques, these models computed the radiative trans- 
fer through spherical shells of dust illuminated by point 
sources. The model used by Gordon et al. (1994) al- 
lowed for multiple illuminating stars arbitrarily distributed 
throughout a spherical nebula and the images of such a 
system could be constructed for arbitrary lincs-of-sight. 
As such, the Gordon ct al. (1994) model is the direct an- 
cestor of the DIRTY model. It is important to note that 
the algorithm used by Gordon et al. (1994), while giv- 
ing the correct total scattered light value, gave too fiat a 
distribution of scattered light across the system. The cor- 
rect algorithm for efficient computation of the appearance 
of a system for a particular line-of-sight is described by 
Yusef-Zadeh, Morris, & White (1984). 

The motivation for the DIRTY model grew out of the 
conference "Dust Survival in Interstellar/Intergalactic Me- 
dia" held at the Space Telescope Science Institute in 1994. 
One of the major concerns at this conference was that all 
Monte Carlo radiative transfer models at that time as- 
sumed a smooth distribution of dust, yet the interstellar 
medium is known to be clumpy over a large range of size 
scales (Colomb, Poppel, & Heiles 1980; Scalo 1990; Rosen 
& Bregman 1995). We constructed the DIRTY model to 
answer this concern. The main change was to move from 
spherical shells to rectangular cells as the basic unit of the 
dust density distribution. Using the simple system of a 
spherical nebula with a central illuminating star, Witt & 
Gordon (1996) explored the differences between the ra- 
diative transfer in smooth and clumpy 3-dimonsional dust 
distributions. The implications of clumpy dust for galax- 
ies were discussed by Witt & Gordon (2000) using the 
concept of spherical galactic environments. The DIRTY 
model has been applied to starburst galaxies (Gordon ct 
al. 1997, 2000), spiral galaxies (Kuchinski et al. 1998), 
the UW Cen reflection nebula (Clayton et al. 1999), and 
the nucleus of M33 (Gordon et al. 1999). 

The DIRTY model was constructed to compute the ra- 
diative transfer of photons from arbitrary distributions of 
photon emitters through arbitrary distributions of dust. 
This model self-consistently models the scattering and ab- 
sorption of photons through dust and the re-emission of 
photons from dust. The algorithm for the radiative trans- 
fer, including polarization, is based on the work of Witt 
(1977), Yusef-Zadeh, Morris, & White (1984), and Code 
& Whitney (1995). The details of this algorithm are the 
subject of this paper. The re-emission of energy absorbed 
by the dust is the subject of a companion paper (Misselt 
et al. 2000a). The radiative transfer of the re-emitted 
photons is handled by the same algorithm as the stellar 



photons. This requires an iterative procedure to account 
for dust absorption and scattering of re-emitted photons 
and their subsequent re-emission. The details of the algo- 
rithms for the dust emission and iterative procedure are 
given by Misselt et al. (2000a). 

2. MONTE CARLO RADIATIVE TRANSFER 

The goal of the DIRTY model is to compute the radia- 
tive transfer of photons through arbitrary distributions of 
dust. The lack of symmetries motivated us to use Monte 
Carlo techniques which are based on the work of Witt 
(1977), Yusef-Zadeh, Morris, & White (1984), and Code 
& Whitney (1995). The forward scattering nature of dust 
grains (Gordon et al. 1994, 1997) and the observed multi- 
phase characteristics of the distribution of dust strongly 
argue for radiative transfer models based on Monte Carlo 
techniques. Such techniques rely on the probabilistic un- 
derstanding of the interaction of a photon with a dust 
grain and sufficient computing power to evaluate this in- 
teraction for many paths through the distribution of the 
dust. This approach allows us to efficiently compute what 
a system of stars and dust will look like for any line-of- 
sight. In addition to the algorithm described in §2.1-2.3, 
the model requires the physical description of the dust 
grain properties, the dust distribution, the distribution of 
photon emitters, and a pseudo-random number generator. 

The physical description of the dust grain properties is 
given by the wavelength dependent behavior of the dust 
optical depth, albedo, and scattering phase function. The 
optical depth determines the probability of a photon in- 
teraction with a dust grain, the albedo determines the 
probability of the dust grain scattering (or absorbing) the 
photon, and the scattering phase function gives the prob- 
ability of the photon being scattered at a specific angle as 
well as the change in the polarization state of the photon 
as a result of the scattering. We have taken the physical 
description of the dust properties from work by Clayton et 
al. (2000) and Kim, Martin, & Hendry (1994). Clayton 
et al. (2000) used MEM techniques to model the size dis- 
tribution of dust grains for dust extinction curves in the 
Milky Way, Large Magellanic Cloud, and Small Magellanic 
Cloud. This work provides albedos and scattering phase 
functions appropriate for the range of known interstellar 
dust extinction curves. 

The DIRTY model allows for arbitrary distributions of 
dust limited only by the amount of memory available to 
store the distribution. This is accomplished by represent- 
ing the dust distribution with a 3-dimensional grid. Each 
grid cell represents a region of uniform dust density. Thus, 
any dust geometry can be represented with the smallest 
scale of inhomogeneity being the size of a grid cell. The 
physical dimensions of each axis of the grid are specified 
by separate 1-dimensional arrays. This allows the phys- 
ical dimensions of one axis to vary making it possible to 
efficiently model both cubical geometries as well as thin 
pancake-like geometries. In addition, this geometry allows 
for the representation of multi-phase or clumpy dust distri- 
butions. In the simplest form, this is a two-phase medium 
with high-density clumps and a low density inter-clump re- 
gion. This type of dust distribution can be described using 
three parameters: the fractional filling factor of the high 
density dust (//), the density ratio between the low and 
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high density dust {k^/ki), and the size of a grid cell in com- 
parison to the system size (l/A^, where N is the number 
of grid cells on an axis). The differences in the radiative 
transfer between homogeneous and two-phase dust distri- 
butions is discussed at length by Witt & Gordon (1996). 

The tracking of the photons through the 3-dimensional 
grid is simple and quite efficient. A photon's location and 
direction of travel is uniquely defined by its x, y, and z 
positions and the u, v, and w direction cosines. When a 
photon enters a grid coll, the side through which it exits 
can be quickly computed in the following manner. Using 
the X, y, and z dimensions of the cell and the direction 
cosines of the photon, the distance the photon would have 
to travel in each of the x, y, and z directions to exit the 
cell through a plane parallel to the wall of the cell in those 
directions can be computed. The photon exits through the 
side which results in the shortest distance traveled. The 
photon's position can be stepped by this distance and an- 
cillary information updated (e.g., optical depth traveled). 

Just as our model permits the use of arbitrary dust dis- 
tributions, it also allows photons to be emitted from ar- 
bitrary source distributions. These can be a single star, 
a distribution of stars, or any definable surface (e.g., an 
accretion disk). The sources can emit isotropically or only 
in specific directions. We use the "ran2" pseudo-random 
number generator which is described by Press, Teukolsky, 
Vetterling, & Flannery (1992) and has a period greater 
than 2 x 10^^. 

2.1. Radiative Transfer Algorithm - A Photon's Life 

The transfer of photons through dust using Monte Carlo 
techniques is historically based on following a single pho- 
ton through a distribution of dust. In this paradigm, the 
chance of the photon being absorbed or scattered is deter- 
mined using a random number generator with the results 
dictating the life of that photon. Therefore, a particular 
photon might exit the system without interacting with a 
dust grain, be absorbed on its first scattering (or subse- 
quent scatterings), or scatter one or more times before 
exiting the system. This is not a particularly efficient 
method of computing the direct, scattered, or absorbed 
light from a system as each photon can contribute to only 
one of these quantities. In addition, each photon exits the 
system in a particular direction which is not necessarily 
the direction from which one wants to observe the system. 
The efficiency of the Monte Carlo method can be increased 
greatly by assigning the photon a weight which allows ev- 
ery photon to contribute to the absorbed energy, the direct 
light, and the scattered light for a particular line-of-sight 
(Witt 1977; Yuscf-Zadeh, Morris, & White 1984). 

The algorithm we use for the radiative transfer in the 
DIRTY model is based on the work of Witt (1977), Yusef- 
Zadeh, Morris, & White (1984), and Code & Whitney 
(1995). In the following description of our algorithm, we 
will use the concept of the life of a photon. To track each 
photon's state, we tabulate its position, direction of travel 
(via direction cosines u, v, and w), weight, and polariza- 
tion (via Stokes I, Q, U, and V values). 

A photon is born according to an input source distri- 
bution and its initial direction is usually chosen from an 
isotropic distribution (e.g., eq. 5-7 of Witt (1977)). The 
assumption of isotropic emission can be modified accord- 



ing to the particular problem being studied. For example, 
a mask with different sized holes at specific locations was 
applied in the DIRTY model to the central star of the 
nebula surrounding the R Corona Borealis star UW Cen 
(Clayton et al. 1999). In this case, only the photons emit- 
ted in the direction of the holes in the mask were allowed 
to move to the next step. The initial weight of the photon 
is usually 

= L«, (1) 

where L" is the luminosity associated with the ath pho- 
ton and a nms from 1 to N. Usually L" = L/N where 
L is the luminosity of the source distribution at the mod- 
eled wavelength and N is the number of photons in the 
model run. The value of can be varied depending 
on the particular problem being studied to optimize the 
model running time (Yusef-Zadch, Morris, & White 1984). 
The photon is initially assumed to be unpolarized (i.e. 
5o" = (Jo«,Q^,[/o",^o") = (l,0,0,0)). 

The flux corresponding to the part of the photon which 
escapes from the system in the direction of the observer is 

^ag-r(o6«)? 1 ^ (2) 



47rd2' 

where T(o6s)g is the optical depth along the path from the 
birth position of the photon to the surface of the system 
in the direction of the observer and d is the distance to 
the system being modeled. In general, r(o6s)o is different 
for each photon emitted imless the system being modeled 
has a single central point source embedded in a spherical, 
homogeneous dust distribution. 

The first scattering of the photon is forced to insure that 
every photon contributes to the scattered light (Cashwell 
& Everett 1959). The optical depth to the first scattering 
site is 



■In 



l-^(l-e--") 



(3) 



where ^ is a random number between and 1 and t" is 
the optical depth to the surface of the nebula in the direc- 
tion the ath photon is traveling. The weight of the photon 
after the first scattering is then 

(4) 



Wf = Wo"a (l - e"^") 



where a is the dust albedo. The fraction of the photon 
which is absorbed at the scattering site is 



A'^ = Wo"{l 



1 



(5) 



The flux corresponding to the fraction of the photon which 
is scattered at the flrst scattering toward the observer is 

F« = w^«e--(°''«)?$(^(o6s)?)^, (6) 

where t(o6s)" is the optical depth from the first scatter- 
ing site along the direction towards the observer, <E>(6') is 
the scattering phase function, and 6{obs)f is the angle be- 
tween the direction the photon was traveling before the 
scattering and the direction towards the observer. This 
angle is easily calculated for the nth scattering using 

COs{9{obs)^) = «_iUo5s + v"_-^Vobs + W^_iWobs) (7) 

where u„-i, Wn-i, and w„-i are the direction cosines of 
the photon before it scatters for the nth time and Uofes, 
Vobs, Wobsi are the direction cosines from the nth scatter- 
ing site towards the observer. This equation is equivalent 
to eq. 17 of Yusef-Zadeh, Morris, & White (1984). 
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L(*) = 



(9) 



The polarization state of the flux corresponding to the 
fraction of the photon which reaches the observer from the 
nth scattering is calculated using 

= Lin - i^)R{e{obsm-zns:-i (8) 
where S" is the new Stokes vector, L{—if) rotates S"_i 
from the reference frame to the scattering plane, R{0) is 
the scattering matrix, and L(7r— 12 ) rotates the new Stokes 
vector from the scattering plane to the reference frame. 
The reference frame for the Stokes vector is set to be the 
z-axis. The rotation matrix is 

10 GO 
cos(2*) sin(2*) 
-sin(25') cos(2^') 
1 

where ^ is the rotation angle. The last two equations are 
eqs. 2-3 of Code & Whitney (1995). 

The scattering matrix, R(0), is the 4x4 Mueller matrix 
and is also known as the phase matrix (Bohren & Huff- 
man 1983; Mishchenko, Hovenier, & Travis 2000). For 
this work, we have used calculations of R{6) by Clayton et 
al. (2000) for size distributions of spherical dust grains. 
As a result, there are only 8 nonzero elements of R{9), 
with Rii = i?22) Ri2 = R211 R33 = Rii-, and R34 = —R43. 
The angles i" and i" are calculated from the knowledge 
of the direction the photon is traveling and 9{obs)" using 
spherical geometry. A useful figure for visualizing the cal- 
culation of if and is Fig. 1 of Code & Whitney (1995). 
The Stokes polarization vector is renormalized after each 
scattering so that /" = 1. As a result, the flux contained 
in the Q polarization component is just Q"F". Similar 
equations give the U and V polarization component fluxes. 

The direction into which the photon is scattered at 
the nth scattering is described by the spherical angles 
9" and cp". The 9^ angle is determined from the the 
dust scattering phase function. We use either a Henyey- 
Greenstein phase function or the full scattering matrices. 
The Henyey-Greenstein phase function is 



Hcos{9),g) 



1 



l3/2 



(10) 



47r [l + g"^ -2g 008(9)]- 

which is eq. 3 of Witt (1977) and where g = {cos 9) is 
the scattering phase function asymmetry. For the Henyey- 
Greenstein phase function, the cosine of the scattering an- 
gle is calculated using 



cos 6'" = 



1 

25 



1 
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1-5 + 25^ 



(11) 



which is eq. 19 of Witt (1977). For the full scattering ma- 
trices, the scattering angle is determined from numerically 
inverting the function 

H^n) = i^-iRiiK) + Ql.Mei), (12) 

where i?ii(6'") and i?22(6'") are elements of the scattering 
matrix and and Q^-i are elements of the photon's 
Stokes vector. Unlike Code & Whitney (1995) who ap- 
proximate ^^{9) as i?ii(0), we use the polarization depen- 
dent scattering phase function. The 0^ angle of scattering 
is determined using 

C=7r(2e-1). (13) 
Given the direction into which the photon is scattered [9^ 
and (f)'^), the new direction cosines can be computed from 



the old direction cosines (e.g., eq. 22 of Witt (1977)). The 
polarization state of the photon after the nth scattering is 
computed using eq. 8 with 9"^ replacing 9{obs)'^. 

The subsequent scatterings are not forced. By forcing 
the first scattering we insure that every photon contributes 
to the scattered flux. By not forcing the subsequent scat- 
terings we insure that the calculation is only done for the 
scatterings which contribute significantly to the scattered 
flux. The optical depth to the nth scattering, , is deter- 
mined using 



■ln$ 



(14) 



which holds for n > 2. If r" is greater than the optical 
depth to the surface in the direction the photon is travel- 
ing, the photon escapes the system. If not, the photon is 
scattered in a direction determined using eqs. 11-13. The 
weight of the photon after the nth scattering is 



(15) 



which holds for n > 2. The fraction of the photon which 
is absorbed at the scattering site is 



AZ = {1- a)W:_ 



(16) 



which holds for n > 2. The fraction of the photon which 
is scattered towards the observer is 



$((^(o6s)^) 



1 



(17) 



which holds for n > 1. The total number of scatterings 
a photon undergoes is limited to be less than Umax- For 
low to moderate optical depths, we set to n„-iax = 50 to 
avoid spending significant computation time on photons 
whose weight is very small as oc a". It is important to 
remember that there are times (e.g., high r values) when 
the value of rimax should be set to a higher value. For 
example, if it is important to know the energy absorbed 
in very high optical depth parts of a dust distribution the 
value of Umax should be increased. 

2.2. Integrated Quantities 

While the previous section has described in detail the 
life of a single photon, the results we are interested in are 
integrated quantities describing the system being modeled. 
These integrated quantities are images of the system (e.g., 
direct and scattered light images) and the 3-dimensional 
matrix of the energy which was absorbed by the dust. 

Images of the system are constructed by projecting the 
birth position (for the direct light image) or the nth scat- 
tering position (for the scattered light image) of the pho- 
ton onto the sky. In the DIRTY model, we embed the dust 
density grid in a sphere making the sky a plane which is 
tangent to the sphere in the direction of the observer. This 
plane is divided into rectangular pixels and the images are 
built up as each photon is run through the model. Special 
care must be applied when constructing the images which 
give the polarization state of the scattered flux, specifl- 
cally the Q and U images. For the radiative transfer of 
the photon, the polarization state has been referenced to 
the z-axis of the system. As part of the construction of 
these images, the polarization vector of the photon must 
be rotated so that it is referenced to the y-axis of the image 
plane. 
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The direct light surface brightness image is constructed 



using 



a 



'pixel 



(18) 



where is the direct flux for the ath photon (see eq. 2), 
^pixei is the surface area of pixel in stcradians, and 
the sum over a is done only for photons whose birth posi- 
tion fall within pixel when projected onto the sky. 

The scattered light surface brightness image is con- 
structed using 

^5(^'^)=jY^EE^" (19) 

^ an 

where is the scattered flux for the nth scattering of 
the ath photon (sec cqs. 6 and 17), and the sums over a 
and n are done only for photons and scatterings of that 
photon whose scattering sites fall within pixel when 
projected onto the sky. The images which describe the po- 
larization state of the scattered light image are constructed 
by first rotating the 5" Stokes vector from the coordinated 
system referenced to the z-axis to a coordinate system ref- 
erenced to the y-axis of the image plane. This is done 
using 

s'^ = Lins: (20) 

where /?" is the angle which rotates the Stokes vector be- 
tween the two coordinate systems. The Q image, and sim- 
ilarly the U and V images, are constructed using 

Q(^,j) = o^EE«n^n (21) 
S ''pixel ^ ^ 

where the sums over a and n are done over the same limits 
as eq. 19. Of course, the sum of the direct and scattered 
light images gives the image one would observe at a tele- 
scope. 

The 3-dimensional absorbed energy matrix is calculated 
using 

a n 

where the sum over a and n are only done over scatterings 
which happen in cell k) and is the energy which 
is absorbed the nth scattering site of the ath photon (sec 
eqs. 5 and 16). The 3D absorbed energy matrix is what 
is needed to compute the dust emission spectrum (Misselt 
et al. 2000a). 

2.3. Uncertainties 

Due to our use of photon weights, the uncertainties in 
output quantities (i.e., scattered intensity, polarization, 
etc.) can not be c;omputcd directly from the square root 
of the number of photons as is usually done when non- 
weighted Monte Carlo techniques are used. The ability to 
calculate uncertainties is crucial when using Monte Carlo 
techniques as the accuracy of model results are dependent 
on the number of photons run. If such uncertainties can 
be calculated during the model run, they can be used to 
dynamically set the number of photons needed in a model 
run to achieve a user input accuracy. We have adopted a 
simple technique for computing the upper limits on the un- 
certainties in the output quantities. If the output quantity 
is X, then 



N 



(23) 



where x'^ is the contribution of the nth scattering of the 
ath photon to X, M is the total number of number pho- 
tons or scatterings in the model run, and x is the average 
contribution each photon or scattering makes to X. In 
the case of the direct flux, the sum over n is dropped and 
M = N . The uncertainty in X is then 



(24) 



where is the standard deviation of x. The value of ax 
is calculated using 



1 



N 



TyEE«--)^ 



(25) 



Again, in the case of the direct flux, the sum over n is 
dropped and M = N. 

The imcertainty ax is only an upper limit on the true 
uncertainty since this quantity measures both the Monte 
Carlo noise associated with running a finite number of pho- 
tons and the intrinsic variation in the quantity. For exam- 
ple, the total scattered flux from a nebula can be computed 
by directly summing every photon's scattered weight. The 
uncertainty, ax, in this quantity has a contribution from 
the intrinsic variation of the scattered flux across the neb- 
ula as well as the Monte Carlo noise. An upper limit which 
would be closer to the true uncertainty can be computed 
by computing the point-by-point uncertainties in an im- 
age of the nebula. The point-by-point uncertainties can 
be calculated using 



(rx{i,jf = 



1 



M{^,J){M{^,J)-l)^^^< ^^^'^'^^^ 

= M(^(^-^T ^''^ 

where M{i,j) is the number of scatterings contributing to 
output quantity in pixel {i,j), the sum over a and n is 
only done for those photons which contribute to the out- 
put quantity in pixel {i,j), and X{i,j) is calculated using 
eq. 18, 19, or 21. In the case of the direct surface bright- 
ness image, the sum over n is dropped and M{i,j) is the 
number of photons contributing to pixel {i,j). 

The uncertainty in the quantity X can then be calcu- 
lated using 



'X 



Y,<^x{i,jf- 

(id) 



(27) 



Calculating the uncertainty with eq. 27 will result in a 
lower value of the uncertainty because the intrinsic varia- 
tion of the quantity X{i,j) across the nebula will not be 
included in the uncertainty calculation. This calculation 
of the uncertainty is still only an upper limit as the un- 
certainty at a specific point in a nebula will still include 
a contribution from photons having intrinsically different 
weights due to having emerged from different depths in 
the nebula. 

We have tested this method of calculating uncertain- 
ties by running the same model multiple times, but with 
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Fig. 1. — The uncertainties in the scattered flux (a) and Q component of the polarized flux (b) are plotted versus the number of photons 
run. At each value of the number of photons, three values of the uncertainties are plotted. The solid line plots the uncertainty as calculated 
from the variation in the 100 difTerent runs, the dashed line plots the average uncertainty as calculated for the total quantity using eq. 25, 
and the dot-dashed line plots the average uncertainty as calculated using eq. 27. The value, averaged over the 100 10® photon runs, of the 
scattered flux (a) is 0.237 and the Q component of the polarized flux (b) is —3.6 x 10~®. These numbers are for the case where L = 1 and 
d = 1. Note that the Q component of the polarized flux is equivalent to zero within the uncertainty, (—3.6 ± 5) x 10~®, as expected for this 
spherically symmetric model. 



a different random number seed each time the model was 
run. This allowed us to compute the true uncertainty in an 
output quantity directly from the variation of the quantity 
between model runs. By comparing the true value of the 
uncertainty with the values computed using eqs. 25 and 
27, we were able to evaluate the accuracy of our simple 
method of estimating uncertainties in output quantities. 
For this testing, wc chose to use a sphere with a homo- 
geneous dust distribution, a ry = 1, V band Milky Way 
dust grain properties, and a central illuminating star (see 
Fig. 2a). Other optical depths give similar results. We 
ran the model 100 times and varied the total number of 
photons between 10^ and 10"''. Figure 1 displays the un- 
certainties in the scattered flux and Q component of the 
polarized flux as a function of the number of photons run. 
The general trend is for the uncertainty calculated using 
eq. 25 or eq. 27 to underestimate the actual uncertainty by 
smaller amounts as the number of photons run increases. 
This is due to small number statistics, especially when the 
uncertainty was calculated using eq. 27. For a large num- 
ber of photons (e.g., 10^), the uncertainty calculated using 
eq. 25 or 27 is a very good estimate of the actual uncer- 
tainty. The reason is that the majority of the scattered 
flux and Q component of the polarized flux comes from 
the central region of the nebula (see Fig. 2a) and the in- 
trinsic variation of the scattered flux in the central region 
is small. This implies that the uncertainty calculated from 
eq. 25 or 27 is dominated by Monte Carlo noise for this 
model. There are model systems where this will not be 
the case and, as a result, care must be taken calculating 
the uncertainty using the method outlined above. 

2.4. Comparison with Other Models 

We tested the results of the DIRTY model against those 
produced by Monte Carlo radiative transfer models which 
do not weight photons and models which use the Witt 
(1977) photon weighting. These models include ones which 
we have coded as well ones others have coded (J. Bjork- 
man 1999, private communication; K. Wood 1999, private 



communication). For computational reasons, these mod- 
els are usually restricted to spherically symmetric systems 
with smoothly varying radial dust distributions. Our two 
main test cases were for r = 1 and r = 10. We adopted an 
albedo of 0.6 and a scattering phase function asymmetry of 
0.6. In all cases, the Yusef-Zadeh, Morris, & White (1984) 
photon weighting method produced statistically similar re- 
sults to the other two weighting methods. In addition, we 
computed the wavelength dependence of the polarization 
for active galactic nucleus models similar to those used by 
Manzini & di Serego Alighieri (1996) and found qualita- 
tive agreement with their results. Quantitative agreement 
is more difflcult to test as we used a different dust grain 
model than Manzini & di Serego Alighieri (1996). 

Figure 2 illustrates the images produced by the DIRTY 
model. Figure 2a shows how a spherical nebula with a cen- 
tral illuminating star would look in the V band assuming 
Milky Way type dust with a homogeneous distribution and 
a radial ry = 1. Figure 2b shows how a biconical nebula 
or active galactic nucleus inclined by an angle of 30° would 
look in the V band assuming Milky Way type dust with a 
homogeneous distribution and a Ty — 1. Figure 3 displays 
the spectral energy distribution before and after the inclu- 
sion of dust in a simple starburst system. This is a good 
illustration of how the dust redistributes energy from the 
ultraviolet to the infrared. In addition, the three compo- 
nents (thermal equilibrium, thermal non-equilibrium, and 
aromatic feature emissions) of the dust emission spectrum 
are shown. Starburst systems are investigated in detail in 
Misselt et al. (2000a). 

3. SUMMARY AND DISCUSSION 

We have presented the DIRTY radiative transfer model 
in this paper and a companion paper (Misselt et al. 
2000a). This model uses Monte Carlo techniques to com- 
pute the radiative transfer of photons emitted from ar- 
bitrary stellar distributions through arbitrary dust dis- 
tributions. The weighting scheme used in the DIRTY 
model results in the eflacient computation of the appear- 
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Fig. 2. — Example images of the DIRTY model are shown for a spherical system (a) and a conieal system (b). The spherieal system is an 
example of a refleetion nebula with an embedded star. The conical system is an example of a bipolar nebula or a quasi-stellar object with 
an opening angle of 45° and an inclination of 30° with the bottom lobe tilted toward the observer. The grayscale displays the distribution of 
scattered light and the line segments the strength and orientation of the polarization. The maximum length of the line segments corresponds 
to a polarization of 48% (a) and 56% (b). 




0.1 1.0 10.0 100.0 

X [yLim] 

Fig. 3. — The spectral energy distribution (SED) of a starburst stellar population is plotted with and without dust. The stellar population 
is for constant star formation of 1.6 MQ/year for a duration of 100 million years. The dust is Milky Way-like with a = 1- The dust 
distribution is the clumpy SHELL geometry (see Misselt et al. (2000a) for additional details). 
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ance of such systems from any viewing direction. The dust 
re-emission inchidcs cquihbrium thermal emission, non- 
equihbrium thermal emission, and the aromatic feature 
emission. The dust emission is computed self-consistently 
with the dust absorption and scattering using an iterative 
technique. The details of the dust emission are described 
by Misselt et al. (2000a). 

The DIRTY model can be used to study a wide range 
of astrophysical objects due to its flexibility. The gen- 
eral behavior of multi-phase dust has been studied us- 
ing the DIRTY model in reflection nebula environments 
(Witt & Gordon 1996) and galactic environments (Witt 
& Gordon 2000) . The very non-symmetric reflection neb- 
ula surrounding the R CrB star UW Gen was successfully 
modeled using the DIRTY model (Clayton et al. 1999). 
The optical depth of spiral galaxies was investigated us- 
ing an early version of the DIRTY model (Kuchinski et al. 
1998). The finding that the lack of a 2175 A depression 
in the spectra of UV-selected starburst galaxies was not 
due to radiative transfer effects but due to the dust grain 
properties was accomplished with the DIRTY model (Gor- 
don et al. 1997). The ultraviolet through near-IR spectral 
energy distribution of the M33 nucleus was modeled using 
the DIRTY model and discovered to be a post-starburst 
stellar population surrounded by a shell of Milky Way-like 
dust (Gordon et al. 1999). The general behavior of the 
infrared emission of starburst galaxies has been studied by 
Misselt et al. (2000a). The full ultraviolet through far- 
infrared spectral energy distributions of a handful of star- 
burst galaxies was studied by Misselt, Gordon, & Clayton 
(2000b). The range of objects which the DIRTY model 
has already been applied to is a good illustration of the 
usefulness of such a radiative transfer model. 



In addition to continuing to use the DIRTY model to 
explore the nature of dust and its affects on a variety of 
astrophysical systems, we plan to continue to improve the 
model itself. One such improvement would be to modify 
the current way the dust distribution is stored to allow 
for grid points to be subdivided independently allowing 
for a larger range of size scale to be efficiently modeled 
(Wolf, Henning, & Stecklum 1999). For example, this 
would allow for modeling the physically small region near 
the accretion disk in a QSO while simultaneously mod- 
eling the large scattering regions in the jets of the QSO. 
Another improvement would be to include electron scat- 
tering which would improve the application of the DIRTY 
model to the study of QSOs. Finally, the code could be 
parallelized to allow larger dust distribution grids to be 
studied in shorter amounts of time. 

As there are always more interesting astrophysical sys- 
tems to be stiidicd than time in the day, anyone interested 
in using the DIRTY model should contact the authors for 
possible collaboration. 
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